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For spiking neural networks we consider the stability problem of global synchrony, arguably the 
simplest non-trivial collective dynamics in such networks. We find that even this simplest dynamical 
problem - local stability of synchrony - is non-trivial to solve and requires novel methods for its 
solution. In particular, the discrete mode of pulsed communication together with the complicated 
connectivity of neural interaction networks requires a non-standard approach. The dynamics in the 
vicinity of the synchronous state is determined by a multitude of linear operators, in contrast to a 
single stability matrix in conventional linear stability theory. This unusual property qualitatively 
depends on network topology and may be neglected for globally coupled homogeneous networks. 
For generic networks, however, the number of operators increases exponentially with the size of the 
O ■ network. 

We present methods to treat this multi-operator problem exactly. First, based on the Gershgorin 
and Perron-Frobenius theorems, we derive bounds on the eigenvalues that provide important infor- 
mation about the synchronization process but are not sufficient to establish the asymptotic stability 
or instability of the synchronous state. We then present a complete analysis of asymptotic stability 
for topologically strongly connected networks using simple graph-theoretical considerations. 
■ For inhibitory interactions between dissipative (leaky) oscillatory neurons the synchronous state 

/"*\ ' is stable, independent of the parameters and the network connectivity. These results indicate that 

pulse-like interactions play a profound role in network dynamical systems, and in particular in 
. the dynamics of biological synchronization, unless the coupling is homogeneous and all-to-all. The 

concepts introduced here are expected to also facilitate the exact analysis of more complicated 
dynamical network states, for instance the irregular balanced activity in cortical neural networks. 
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I. INTRODUCTION 

(N 
t> 

t^J- ■ Understanding the collective dynamics of biological neural networks represents one of the most challenging problems 
in current theoretical biology. For two distinct reasons, the study of synchronization in large neuronal networks plays 
a paradigmatic role in theoretical studies of neuronal dynamics. Firstly, synchrony is an ubiquitous collective feature 
of neural activity. Large-scale synchronous activity has been observed by spatially coarse methods such as electro- 
or magnetoencephalography. Complementing experiments of parallel recordings of spiking activity of individual cells 
have shown that the synchronization of firing times of indiviual units is often precise with a temporal scatter of the 
order of a few milliseconds [HJ3|. This precise locking has been observed over significant distances in the cortex 
and even across hemispheres [3|. Synchronous activity also plays an important role in pathological state such as 
epileptic seizures Secondly, the synchronous state is arguably the simplest non-trivial collectively coordinated 
state of a network dynamical system. Mathematically, it is therefore one of the most throroughly investigated states 
in the dynamics of biological neural networks. Following the paradigm set by the seminal works of Winfree and 
Kuramoto on the dynamics of biological oscillators, most studies of synchronization have utilized either temporal or 
population averaging techniques to map the pulse-coupled dynamics of biological neural networks to effective models 
of phase-coupled oscillators or density dynamics @, 0,0, B B, ES EI El El EI El, El, E3, El, El, US HH, SI . 

While this approach has proven to be very informative, it provides exact results only under highly restrictive 
conditions. In facj^studies approaching the dynamics of biological neural networks using exact methods 0, US US 
US US US US US US El, HI E3, HH, have over the past decade revealed numerous examples of collective behaviours 
that obviouslly are outside the standard repertoir of behaviours expected from a smoothly coupled dynamical system. 
These include several new and unexpected phenomena such as unstable attractors (23l . [HI, EH] , stable chaos H 22, 
l36| . topological speed limits to coordinating spike times [U, and extreme sensitivity to network topology [33]. 
The occurence of these phenomena may signal that the proper theoretical analysis of collective neuronal dynamics 
mathematically represents a much more challenging task than is currently appreciated. 

In this article, we study the impact of pulse-coupling, delayed interactions and complicated network connectivity 
on the exact microscopic dynamics of neural networks and expose the mathematical complexity that emerges already 
when considering the seemingly simple problem of neuronal synchronization. Utilizing the Mirollo-Strogatz phase 
representation of individual units, we present an analytical treatment of finite networks of arbitrary connectivity. 
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The results obtained unearth an unanticipated subtlety of the nature of this stability problem: It turns out that 
a single linear operator is not sufficient to represent the local dynamics in these systems. Instead, a large number 
of linear operators, depending on rank order of the perturbation vector, is needed to represent the dynamics of 
small perturbations. We present methods to characterize the eigenvalues of all operators arising for a given network. 
Universal properties of the stability operators lead to exact bounds on the eigenvalues of all operators and also provide 
a simple way to prove plain but not asymptotic stability for any network. We then show for topologically strongly 
connected networks that asymptotic stability of the synchronized state can be demonstrated by graph-theoretical 
considerations. We find that for inhibitory interactions the synchronous state is stable, independent of the parameters 
and the network connectivity. A part of this work that considers plain (non-asymptotic) stability without networks 
constraints has been briefly reported before for the case of inhibitory coupling [3a]. These results indicate that pulse- 
like interactions play a profound role in network dynamical systems, and in particular in the dynamics of biological 
synchronization, unless the coupling is homogeneous and all-to-all. They highlight the need for exact mathematical 
tools that can handle the dynamic complexity of neuronal systems at the microscopic level. 



II. THE PHASE REPRESENTATION OF PULSE-COUPLED NETWORKS 



We consider networks of N pulse-coupled oscillatory units, neurons, with delayed interactions. A phase-like variable 
4>j(t) 6 (— oo, 1] specifies the state of each oscillator j at time t such that the difference between the phases of two 
oscillators quantifies their degree of synchrony, with identical phases for completely synchronous oscillators. The free 
dynamics of oscillator j is given by 



Whenever oscillator j reaches a threshold 



the phase is reset to zero 



M*) = 1 (2) 



= (3) 



and a pulse is sent to all other post-synaptic oscillators i <E Post(j). These oscillators i receive this signal after a delay 
time t. The interactions are mediated by a function U (4>) specifying a 'potential' of an oscillator at phase <j>. The 
function U is assumed twice continuously differentiable, monotonically increasing, 

U' > 0, (4) 

concave (down), 

U" < 0, (5) 

and normalized such that 

17(0) = and U{1) = 1. (6) 

For a general U {4>) we define the transfer function 

H E (<j>)=U- 1 (U(c(>)+e) (7) 

that represents the response of an oscillator at phase 4> to an incoming sub-threshold pulse of strength e. Depending 
on whether the input received by oscillator % from j is sub-threshold, 

U{4>) + en < 1, (8) 

or supra-threshold, 

[/(0)+ £ii >1, (9) 
the pulse sent at time t (Eq. (JSJ)) induces a phase jump after a delay time r at time t + r according to 

6 .(( f + T ) + )- J H eij (<f>i(t + T)) if tf(&(t + 7-))+ eij -<l ( 

( P^ t + T > >-\0 if U(Mt + r))+e l3 >l ■ [W) 
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Figure 1: An incoming pulse of strength e induces a phase jump (f>f := 4>i(t + ) = U~ 1 {U(4>\{t)) + e) = H e (4>i) that depends 
on the state (f>\ := <j>i(t) of the oscillator at time t of pulse reception. Due to the curvature of U, an excitatory pulse (e > 0) 
induces an advancing phase jump. If the incoming pulse puts the potential above threshold {U{4>\) + e > 1), the phase is reset 
to zero (<f)f = 0). An inhibitory pulse (e < 0) would induce a regressing phase jump such that the phase may assume negative 
values (not shown). 



In the second case also a pulse is emitted by oscillator i . The phase jump (Fig.Q]) depends on the phase 0i(i+r) of the 
receiving oscillator i at time t + r after the signal by oscillator j has been sent at time t, the effective coupling strength 
Eij, and the nonlinear potential function U. The interaction from unit j to unit i is either excitatory (ey > 0) inducing 
advancing phase jumps (cf. Fig. Q]) or inhibitory (ey < 0) inducing retarding phase jumps. If there is no interaction 
from j to i, we have Eij — and . There are two immediate differences between inhibitory and excitatory inputs. 
First, while inhibitory input < is always sub-threshold, excitatory input Eij > may also be supra-threshold 
and thus induce an instantaneous reset to zero according to Eq. lfT0|). Second, in response to the reception of an 
inhibitory pulse, the phases of the oscillators may assume negative values whereas for excitatory coupling they are 
confined to the interval [0, 1]. We remark here that we do not consider the dynamics (and perturbations) of variables 
that encode the spikes in transit, i.e. spikes that have been sent but not yet received at any instant of time. These 
additional variables make the system formally higher-dimensional. However, under the conditions considered below, 
in particular for networks of identical neurons with inhibitory interactions, earlier rigorous work [35j shows that these 
spike time variables lock to the phase variables once all spikes present in the system initially have arrived and every 
neuron has emitted at least one spike thereafter (which takes finite total time). Thus, for our purposes, we consider 
the dynamics and perturbations of the phase variables only. 

Choosing appropriate functions U the dynamics of a wide variety of pulse coupled systems can be represented. In 
particular, any differential equation for an individual neural oscillator of the form 

V i = f(V i ) + S i (t) (11) 

together with the threshold firing condition 

Vi(t) = 1 => V t (t+)=0 (12) 

where 

S i {t)=^e ij 8{t-t j . m ) (13) 

j,m 

is the synaptic current from the network and tj >m is the time of the m th firing of neuron j. As long as the free 
(Si(t) = 0) dynamics has a periodic solution V(t) with period T and negative curvature, the function U can be taken 
as the scaled solution, 

U(<f>) := V(#T). (14) 

Thus a general class of pulse-coupled oscillator networks, defined by Eq. l|lHll3p . can be mapped onto the normalized 
phase representation (Eqs. (fTj)- (fT0l0 . 

In this paper, we consider a system of TV coupled Mirollo-Strogatz oscillators which interact on directed graphs 
by sending and receiving pulses. The structure of this graph is specified by the sets Pre(i) of presynaptic oscillators 
that send pulses to oscillator i. For simplicity we assume no self-interactions, i ^ Pre(i). The sets Pre(i) completely 
determine the topology of the network, including the sets Post(i) of postsynaptic oscillators that receive pulses from 
i. The coupling strength between oscillator j and oscillator i is given by ey such that 



^ if j £ Pre(i) 
Eij — otherwise. 



(15) 
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Thus a connection is considered to be present if the connection strength is non-zero. All analytical results presented are 
derived for the general class of interaction functions U introduced above. The structure of the network is completely 
arbitrary except for the restriction that every oscillator has at least one presynaptic oscillator. 



III. REGULAR AND IRREGULAR DYNAMICS IN PULSE-COUPLED NETWORKS 



The synchronous state in which 



<j>i(t) = <j>o(t) for all i 



(16) 



is arguably one of the simplest ordered states a network of pulse-coupled oscillators may assume. This synchronous 
state exists if and only if the coupling strengths are normalized such that 



jePre(i) 



: e. 



(17) 



One should keep in mind that this state whether stable or not is typically not a global attractor in complex networks. 
To illustrate this point let us briefly consider a specific example. Figures (Fig. [2] and [3| show numerical results 
from a in a randomly connected network of integrate-and-fire oscillators where U((j>) = Uiy(d>) = 1(1 — e~^ TlF ) and 
Tif = ln(I/(7 — 1)) with strong interactions. This network exhibits a balanced state (cf. 119. [39 . 140;]) characterized 
by irregular dynamics. In this balanced state, found originally in binary neural networks [39j, |4Qj, inhibitory and 
excitatory inputs cancel each other on average but fluctuations lead to a variability of the membrane potential and 
a high irregularity in firing times (see also Figures EH,b display sample trajectories of the potentials U(<j>i) 

of three oscillators for the same random network, making obvious the two distinct kinds of coexisting dynamics. 
Whereas in the synchronous state all oscillators display identical periodic dynamics, in the balanced state oscillators 
fire irregularly, asynchronous, and in addition differ in their firing rates. 
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Figure 2: Coexistence of (a) synchronous and (b) irregular dynamics in a random network (N — 400, p = 0.2, / = 4.0, 
e = —16.0, r = 0.035), cf. [38l |. (a),(b): Trajectories of the potential U((j>i) of three oscillators (angular bars: time scale 
(horizontal) At = 8; potential scales (vertical) (a) AU = 8, (b) AU = 2 ; spikes of height AU = 1 added at firing times). 
(c),(d): Distributions (c) p„ of rates and (d) pcv of the coefficient of variation, displayed for the irregular (dark gray) and 
synchronous (light gray) dynamics. Figure modified from [3H |. 



The latter dynamical difference is quantified by a histogram p v of oscillator rates (Fig. Eh) 

-l 



(18) 



the reciprocal values of the time averaged inter-spike intervals. Here the ti >n are the times when oscillator i fires 
the n th time. The temporal irregularity of the firing-sequence of single oscillators i is measured by the coefficient of 
variation 



CV, 



(19) 
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defined as the ratio of the standard deviation of the inter-spike intervals and their average. A histogram pcv of the 
CVi shows that the irregular state exhibits coefficients of variation near one, the coefficient of variation of a Poisson 
process. Such irregular states occur robustly when changing parameters and network topology. 

Figure [3] illustrates the bistability of these qualitatively different states by switching the network dynamics form 
one state to the other by external perturbations. A simple mechanism to synchronize oscillators that are in a state of 
irregular firing is the delivery of two sufficiently strong external excitatory (phase-advancing) pulses that are separated 
by a time At 6 (r, 1), cf. Fig. [31 The first pulse then leads to a synchronization of phases due to simultaneous supra- 
threshold input, cf. ([7]). If there are traveling signals that have been sent but not received at the time of the first 
pulse, a second pulse after a time At > r is needed that synchronizes the phases after all internal signals have been 
received. In this network the synchronous state is not affected by small random perturbations, whereas large random 
perturbations lead back to irregular dynamics (Fig. [3j . These features clearly suggests that the synchronized states 
is an asymptotically stable local attractor. However, with the exception of special cases no exact treatment of this 
proposition exists. Below we will expose that a general treatment of the stability of apparently simple synchronous 
state reveals an unexpectedly complex mathematical setting. 
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Figure 3: Switching between irregular and synchronous dynamics (TV = 400, p = 0.2, / = 4.0, e = —16.0, r = 0.14). Firing 
times of five oscillators are shown in a time window At = 240. Vertical dashed lines mark external perturbations: (i) large 
excitatory pulses lead to synchronous state, (ii) a small random perturbation (|A0<| < 0.18) is restored, (iii) a sufficiently large 
random perturbation (|A<^| < 0.36) leads to an irregular state. Bottom: Time evolution of the spread of the spike times after 
perturbation (ii), total length At — 0.25 each. Figure modified from [38]. 



IV. CONSTRUCTING STROBOSCOPIC MAPS 



We perform a stability analysis of the synchronous state 

4>i{t) = <po(t) for all i 



(20) 



in which all oscillators display identical phase-dynamics <po (t) on a periodic orbit such that 0o (t + T) — 0o (t) . The 
period of the synchronous state is given by 



where 



T=T+l-a 



a = U -1 (E/"(t) +e). 



(21) 



(22) 



For simplicity, we consider the cases where the total input e is sub-threshold, U(t) + e < 1 such that a < 1. A 
perturbation 



8(0)=:S = (Sx,...,S N ) 

to the phases is defined by 

Si = fo(p) - <h(p) ■ 

If we assume that the perturbation is small in the sense that 

max Si — min Si < r 



(23) 
(24) 
(25) 
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this perturbation can be considered to affect the phases of the oscillators at some time just after all signals have been 
received, i.e. after a time t > t + r if all oscillators have fired at t = t . This implies that in every cycle, first each 
neuron sends a spike before any neuron receives any spike. Such a perturbation will affect the time of the next firing 
events because the larger the perturbed phase of an oscillator is, the earlier this oscillator reaches threshold and sends 
a pulse. In principle, there one may consider other perturbations, in which, for instance, a pulse is added or removed 
at a certain time. As such perturbations are not relevant for questions of asymptotic stability, we do not consider 
them here. 

To construct the stroboscopic period- T map, it is convenient to rank order the elements Si of S in the following 
manner: For each oscillator i we label the perturbations Sj of its presynaptic oscillators j G Pre(i) (for which Eij ^ 0) 
according to their sizes 

Ai,i > A u2 > . . . > A iM (26) 

where 

ki := |Pre(t)| (27) 

is the number of its presynaptic oscillators, called in-degree in graph theory. The index n 6 {1, . . . , ki} counts the 
received pulses in the order of their successive arrival. Thus, if j n = j n (i) & Pre(i) labels the presynaptic oscillator 
from which i receives its n th signal during the period considered, we have 

Aj )Tl = <5j„(i) ■ (28) 

For later convenience, we also define 

A,. = 6, . (29) 

For illustration, let us consider an oscillator % that has exactly two presynaptic oscillators j and j' such that 
Pre(i) = {j,f} and ki = 2 (Fig. Hk,d). For certain perturbations, oscillator i first receives a pulse from oscillator f 




Figure 4: Two signals arriving almost simultaneously induce different phase changes, depending on their rank order. The figure 
illustrates a simple case where Pre(i) = {j, j'} and 5, = 0, (a)-(c) for 5ji > Sj and (d)-(f) for 8j > 5ji . (a), (d) Local patch 
of the network displaying the reception times of signals that are received by oscillator i. Whereas in (a) the signal from j' 
arrives before the signal of j, the situation in (d) is reversed, (b), (e) Identical coupling strengths induce identical jumps of the 
potential U but (c),(f) the phase jumps these signals induce are different and depend on the order of the incoming signals. For 
small \5i\ 1, individual phase jumps are encoded by the Pi >n - 

and later from oscillator j. This determines the rank order, Sj> > Sj, and hence A^i = Sj> and A^2 = Sj (Fig. H£i). 
Perturbations with the opposite rank order, Sj > Sj', lead to the opposite labeling, A^i = Sj and A^2 = Sj' (Fig.[4]i). 

We now consider a fixed, arbitrary, perturbation, the rank order of which determines the A,- in according to the 
inequalities J26|). Using the phase shift function H E {4>) (see J7])) and denoting 



(30) 
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for n G {1, . . . , fcj} we calculate the time evolution of phase-perturbations Si satisfying J25|). Without loss of generality, 
we choose an initial condition near (fo(0) = r/2 . The stroboscopic time-T map of the perturbations, Si h- > <5,(T), is 
obtained from the scheme given in Table HI The time to threshold of oscillator i, which is given in the lower left of 



t 







§ + Si =: | + Ai, 




(7' 1 (t/(r + A,i)+e lJ1 )=:ft,i 


i-A,, 2 


[/-^{/(A.i + +£ii 2 ) =: Pi,2 




U- l (U(j3i tki -i + D iM ) +e ijhi ) =: 




reset: 1 i— > 



Table I: Time evolution of oscillator i in response to ki successively incoming signals from its presynaptic oscillators j n , 
n £ {1, ...,ki}, from which i receives the n th signal during this period. The right column gives the phases <fii(i) at times t 
given in the left column. The time evolution is shown for a part of one period ranging from cj>i w r/2 to reset, 1 — > 0, such that 
0i = in the last row. The first row gives the initial condition 4>i(0) = r/2 + Si . The following rows describe the reception 
of the ki signals during this period whereby the phases are mapped to /3;, n after the 71 th Si gnal has been received. The last 
row describes the reset at threshold such that the respective time T^ ' = r/2 — A;.^ + 1 — Pi,hi gives the time to threshold of 
oscillator i. 



the scheme, 

?f ) -=\- A*, fel + 1 ~ Pi, ki (31) 
is about </>o(0) = r/2 smaller than the period T. Hence the period-T map of the perturbation can be expressed as 

Si (T) = T - - ~ = p iM - a + A lM (32) 

where a is given by Equation lf22j) . 

Equation f32|) defines a map valid for one particular rank order of perturbations. In general, the perturbations of 
all ki presynaptic oscillators of oscillator i lead to fcj! different possibilities of ordering. Thus the number of possible 
maps, fi, is bounded by 

fmaxfci)! </x< (iV-i)!. (33) 

Here the minimum is assumed, [i = (max^ fcj)!, if only one oscillator has exactly fi presynaptic oscillators and all other 
oscillators have exactly one presynaptic oscillator. The maximum, fi — (N — 1)! , is assumed if the oscillators are 
coupled all-to-all such that all connections are present, £y ^ for all i and j ^ i. If the coupling is all-to-all and 
in addition homogeneous, — e/(N — 1) for all i and j / i, all maps are equivalent in the sense that for any pair 
of maps there is a permutation of oscillator indices that transforms one map onto the other. For general network 
connectivities, however, there is no such permutation equivalence. For instance, in random networks with N vertices 
and edges that are independently chosen with identical probability p, the number of maps increases strongly with N, 



V. MULTI-OPERATOR DYNAMICS OF SMALL PERTURBATIONS 



In order to perform a local stability analysis, we consider the first order approximations of the maps derived in the 
previous section. Expanding /J^fc, for small Z\ n <C 1 one can proof by induction (see Appendix) that to first order 

ki 

0i,ki = a + /,Pi,n-lDi,n (34) 
n=l 

where 



U'{U- x {U{t)+e)) 



(35) 
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for n 6 {0, 1, . . . , ki} encodes the effect of an individual incoming signal of strength £,j n . The statement x = y means 
that x — y + J2i n 0(D? n ) as all D iyJl — > 0. Substituting the first order approximation (|34|) into (|32|) using (f30|l leads 
to 



such that after rewriting 



&ii T ) = / ]pi,n-l(Ai,n-l ~ A;, n ) + A 4>fci (36) 



£j(T) = Pi,oAj,o + ^2(Pi,n -Pt,n-l)Ai,„ (37) 




to first order in all Aj in . Since Aj ; „ = flj n (i) for n e {1, . . . , fc^} and A^o = <5j according to Eqs. |28|) and i[29|) . this 
represents a linear map 

S(T) = AS (38) 

where the elements of the matrix A are given by 

- pi,„_i if j = j„ E Pre(i) 

if j = * (39) 
if j i Pre(i) U {i}- 

Because j„ in ([39| identifies the n th pulse received during this period by oscillator i, the first order operator depends 
on the rank order of the perturbations, A = A(rank(J)). The map AS consists of a number of linear operators, the 
domains of which depend on the rank order of the specific perturbation. Thus is piecewise linear in S. This map 
is continuous but not in general differentiable at the domain boundaries where Si — Sj for at least one pair i and j 
of oscillators. In general, signals received at similar times by the same oscillator induce different phase changes: For 
the above example of an oscillator i with exactly two presynaptic oscillators j and j' and equal coupling strengths, 
— , the first of the two received signals has a larger effect than the second, by virtue of the concavity of U ((j)). 
For small |<5j| -C 1, this effect is encoded by the pi <n (see Eq. (|35| and the Appendix for details). Since the matrix 
elements Eq. I|39p are differences of these Pi >n the respective matrix elements Aij and Aij> have in general different 
values depending on which signal is received first. This is induced by the structure of the network in conjunction 
with the pulse coupling. For networks with homogeneous, global coupling different matrices A can be identified by 
an appropriate permutation of the oscillator indices. In general, however, this is impossible. 

VI. BOUNDS ON THE EIGENVALUES 

The above consideration establish that for many network structures when considering the stability of the syn- 
chronous state one is faced with the task of characterizing a large number of operators instead of a single stability 
matrix. Fortunately it is possible to characterize spectral properties common to all operators by studying bounds on 
their eigenvalues. 

A. General properties of matrix elements 

It is important to observe that the matrix elements defined by Eqs. lf38j) and (|39"1) are normalized row-wise, 

AT N 

= Aa + J2 A ij (40) 

ki 



An + A ijn (41) 



n=l 

ki 



Pi,0 + / l(pi,n - Pi,n-l) (42) 

71=1 

Pi,ki (43) 

1 (44) 
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for all i. Here the second last equality holds because the telescope sum equals ViM — Pi.o ■ Therefore, every matrix A 
has the trivial eigenvalue 

Ai = 1 (45) 

corresponding to the eigenvector 

vi = (1,1,.-., 1) T (46) 
reflecting the time-translation invariance of the system. In addition, the diagonal elements 

A ^^ = Ui{U-^ul) + e)) =-- A » (47) 

are identical for all i. Since U is monotonically increasing, U'(cj)) > for all (p, the diagonal elements are positive, 

A > . (48) 

One should note that the matrices A have the properties l(40|) - ()48| independent of single neuron parameters, the 
network connectivity, and the specific perturbation considered. Due to the above properties of the stability matrices, 
bounds on the eigenvalues of a specific matrix A can be obtained from the Gershgorin theorem (see also [H]). 

Theorem VI. 1 (Gershgorin) Given an N x N matrix A — (j4y) and disks 

N 

Ki-.= {zeC\ \z-Au\ <J2\ A a\} ( 49 ) 

for i G {!,..., N}. Then the union 



N 

K := |J Ki (50) 

i=l 

contains all eigenvalues of A. 
Let us remark that for real matrices A the disks Ki in the complex plane are centered on the real axis at An = Aq. 

B. Eigenvalues for inhibitory coupling 

As an application of this theorem to the above eigenvalue problem, we consider the class of networks of inhibitorily 
coupled oscillators, where all < and e < 0. In these cases, all nonzero matrix elements Aij are positive: Since 
U(<fi) is monotonically increasing, U' > 0, and concave (down), U" < 0, its derivative U' is positive and monotonically 
decreasing. Thus all Pi >n (Eq. j[35|)) are positive, bounded above by one, 

< Pi,n < 1, (51) 

and increase with n, 

Pi,n-i < Pi.n ■ (52) 
Hence the nonzero off-diagonal elements are positive, Ay„ = Pi tU — Pi, n -i > such that 

A i:j > (53) 

for alH, j S {1, . . . , N} and 

Q < A < 1. (54) 
As a consequence, for inhibitorily coupled oscillators, 

N N 

^2\A tj \= Y / A ij -A ii = l-A (55) 
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Figure 5: Eigenvalues in the complex plane for inhibitory coupling. The Gershgorin disk, that contacts the unit disk from 
the inside, contains all eigenvalues of the stability matrices A. Black dots show eigenvalues of a specific stability matrix for 
a network of N = 16 oscillators (in which every connection is present with probability p = 0.25) with coupling strengths 
Eij = e/ki , e = —0.2, r = 0.15, and a particular rank order of the perturbation. 

such that all Gershgorin disks Ki are identical and the disk 

Ki = K = {z G C| \z-A \ < 1 -A Q } (56) 

contains all eigenvalues of A. This disk K contacts the unit disk from the inside at the trivial eigenvalue z = \\ = 1 
(Fig. El. 

Since the unit circle separates stable from unstable eigenvalues and structural perturbations may move any but the 
trivial unit eigenvalue, it is interesting to examine whether the unite eigenvalues may be degenerate. For strongly 
connected networks the Perron-Frobenius theorem implies that this eigenvalue is unique demonstrating the structural 
stability of the confinement of eigenvalues to the unit circle. 

If the network is strongly connected, the resulting stability matrix A is irreducible such that the Perron-Frobenius 
theorem [H, 0, [H, Hi| (see also [13, 0|) is applicable (we only state the theorem partially). 

Theorem VI. 2 (Perron-Frobenius) Let A be an N x N irreducible matrix with all its elements real and non- 
negative. Then 

• A has a real positive eigenvalue X max , the maximal eigenvalue, which is simple and such that all eigenvalues Aj 
of A satisfy \Xi\ < X max - 

The Perron-Frobenius theorem implies, that the eigenvalue that is largest in absolute value, here Ai = 1, is unique for 
strongly connected networks. The Gershgorin theorem guarantees that eigenvalues A of modulus one are degenerate, 
A = Ax = 1. Taken together, for strongly connected networks, the non-trivial eigenvalues Ai satisfy 

I Ai| < 1 (57) 

for i e {2, . . . , N}. This suggests that the synchronous state is stable for inhibitory couplings. As pointed out below 
( Sec. IVIIl a proof of stability, however, requires further analysis. 

C. Eigenvalues for excitatory coupling 

Let us briefly discuss the case of excitatorily coupled oscillators, where all > and e > 0. Here the analysis 
proceeds similar to the case of inhibitory coupling. Due to the monotonicity and concavity of U ((/)), we obtain 

P l ,n > 1 (58) 

as well as a decrease with n, 



Pi,n—1 ^ Pi,n 7 



(59) 



such that Aij n = pi >n — Pi. n -i < and thus 

for i <E {1, . . . , N} and j ^ i and 

Consequently, for excitatorily coupled oscillators, 

JV 



Ay < 



A u = A > 1. 
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(60) 
(61) 

(62) 



3=1 



3=1 



such that again all Gershgorin disks K\ are identical and 

Ki = K = {z£C\ \z-Aq\ <A -1}. (63) 
Since Aq > 1, the disk if contacts the unit disk from the outside at z = Ai = 1 (Fig. [6]). If the network is strongly 

-10 12 




Re A 

Figure 6: Eigenvalues in the complex plane for excitatory coupling. Gershgorin disk, that contacts the unit disk from the 
outside, contains all eigenvalues of the stability matrices A. Black dots show eigenvalues of a specific stability matrix for 
a network of N = 16 oscillators (in which every connection is present with probability p = 0.25) with coupling strength 
Eij = e/ki ,e = 0.2, r = 0.15) and a particular rank order of the perturbation. 

connected, Ai = 1 is again unique by the Perron-Frobenius theorem for irreducible matrices because it is the largest 
eigenvalue (in absolute value) of the inverse A -1 of A, if the inverse exists. This result indicates that the fully 
synchronous state is unstable for excitatory couplings. 



VII. STABILITY 

In Sec. EH we found analytical bounds on the eigenvalues of the stability matrices. However, even if only eigenvalues 
A with A | < 1 are present, a growth of perturbations might seem to be possible because (i) the eigenspaces of 
the (asymmetric) matrices A cannot be guaranteed to be orthogonal, and (ii) in general, different matrices A are 
successively applied to a given perturbation. Due to the non-orthogonality (i) the length \\S\\ of a given perturbation 
vector S might increase during one period. Since the stability matrix and the set of eigenvectors may change due to 
(ii), the length of the perturbation vector might increase in the subsequent period as well. Since this procedure may 
be iterated, the eigenvalues Ai of the stability matrices, although satisfying |Af| < 1 for i G {2, ... , JV} and Ai = 1, 
are not sufficient to ensure the stability of the synchronous state. 

Thus, the eigenvalues of the dynamically changing stability matrices guide the intuition about the stability of the 
synchronous state as well as about the speed of convergence (in case of stability) or divergence (in case of instability) . 
Nevertheless, stability cannot be directly inferred from the set of eigenvalues. In the following we illustrate the final 



proof of stability in the simple case of inhibitory coupling where all Ay > ([53 
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A. Plain stability for inhibitory coupling 

For inhibitory networks, the proof of plain (non-asymptotic) stability is simple. Given the fact that for inhibition 
SjLi Aij = 1 an d > 0, the synchronous state is stable because a given perturbation 6 = 6(0) satisfies 

||<S(T)|| := max\Si(T)\ (64) 



/v 



In this section, we use the vector norm 



< max^lAyll^l (66) 

3 

< max 7 \Aij\ max |£J (67) 

i 

= max A$ j max 5^ | (68) 

i 

= max|i5fc| (69) 

k 

= \\S\\. (70) 
:= max J*. (71) 



Thus the length of a perturbation vector can not increase during one period implying that it does not increase asymp- 
totically. This result is independent of the connectivity structure of the network, the special choice of parameters, 
£ij < 0, t > 0, the potential function U(<j>), and the rank order of the perturbation. 

B. Asymptotic stability in strongly connected networks 

Asymptotic stability can be established for networks of inhibitorily coupled oscillators assuming that the network 
satisfies the condition of strong connectivity. A directed graph is called strongly connected if that every vertex can be 
reached from every other vertex by following the directed connections. Thus a network is strongly connected, if every 
oscillator can communicate with each oscillator in the network at least indirectly. It is clear that in a disconnected 
network only the connected components may synchronize completely in the long-term, but these components can not 
be synchronized by mutual interactions. In the proof given below, we do not consider networks that are disconnected. 
We do also not consider networks that are weakly connected (but not strongly connected) such as two globally coupled 
sub-networks which are linked by unidirectional connections from one sub-network to the other. 

The synchronous state may be characterized by a perturbation 6 = 6(0) that represents a uniform phase shift, 

<5 = civi =ci(l,l,...,l) T (72) 

where c\ el, |ci| < 1, and vi is the eigenvector of A corresponding to the eigenvalue Ai = 1 f46|) . Such a perturbation 
satisfies 

6(T) = A6 = 6 (73) 
for all matrices A — A(rank(5)) independent of the rank order rank(<$) and thus 

|]<5(T)|| = ||«|| . (74) 
Now consider a 6 that does not represent the synchronous state, 

6 + civi (75) 

for all ci £ R. Then one might guess that 

||i(T)||<||«||. (76) 
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We assume that (i) 8 is not the synchronous state f75|) and that (ii) all single-oscillator perturbations are non-negative, 
Si > 0. The latter assumption is made without loss of generality, because otherwise a vector proportional to vi can 
be added such that Si > is satisfied. The assumption (ii) implies that the components of the perturbation stay 
non-negative for all times, because 5i(T) = AijSj > for all i such that 5i(lT) > for all i and all Z G N. 

For convenience, we define the largest component of the perturbation 

S M := max Si (77) 

i 

and the second largest component 

S m := max{(5i | Si < S M } (78) 
such that Sm < Sm ■ We also define the index set of maximal components 

M:={je{l,...,]V}|^4f} (79) 
which is always non-empty. We write j £ M if j G {1, . . . , N} \ M. With these definitions we find 

N 

Si(T) = Y. A n 6 i ( 8 °) 

= X m.- • Y. a -Aj ( 81 ) 

jeu jgM 

AijSM + A V S ™ (82) 
A l jS M - Y, A-ijSm (83) 



< 



■E 

N 



= S M Y A H ~ ~ 6 m) A » ( 84 ) 

3=1 j$M 

= S M - (Sm - 5 m ) ^ A v ■ ( 85 ) 

Hence if J2jgM A-ij > the norm of the perturbation vector decreases in one period, 

||5(T)|| = maxtfi(T) < S M = maxfc = ||«|| . (86) 

i i 

There are, however, also perturbations that imply A y - = for (at least) one specific i and all j ^ M such that 
J2js£M Ai 3 = This is the case if and only if there is an oscillator i that receives input only from oscillators j with 
maximal components, Sj — Sm, and itself has maximal component, Si = Sm ■ So suppose that 

3* G {1, .. .,N} Vj G Pre(i) U {i} : Sj = 5 M , (87) 

i.e. Pre(i) U {i} C M. Then 8i(T) = Sm for this oscillator i and hence 

max^(T) = S M (88) 

i 

such that the norm of the perturbation vector does not decrease within one period. 

Nevertheless, if the network is strongly connected, the norm of the perturbation vector is reduced in at most N — 1 
periods: We know that if and only if there is no oscillator i satisfying l(87|). the norm will be reduced l|86p because 
SjgM Aij > for all i by f85|) . So to have maxi<5i(T) = Sm one needs one oscillator i that satisfies l(87j) . i.e. i and 
Pre(i) have to have maximal components. If the vector norm stays maximal another period, 

max5i(2T) = S M , (89) 

i 

not only all j G {i} U Pre(i) but also all j G Pre(Pre(i)) have to have maximal components, Sj = Sm- Iterating this Z 
times, leads to the condition 

maxjtfj (ZT) = Sm 

(90) 

3i Vj G {i}UPre(i)UPre (2) (i)U...UPre ( ' ) (i) : Sj = 5 M 
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where 

Pre ( ' ) (i) := Pre o Pre o . . . o Pre (i) (91) 

l times 

is the set of oscillators, that is connected to oscillator i via a sequence of exactly I directed connections. 

Since for a strongly connected network the union of all presynaptic oscillators and their respective presynaptic 
oscillators is the set of all oscillators 

{i} U Pre(i) U Pre (2) (i) U . . . U Pre (0 (i) = {1, . . . , N} (92) 

for I > N — 1, this leads to the conclusion that 

max ( 5 l ((7V - 1)T) = S M => Vj : & = S M (93) 

such that 



S = S M V! (94) 

contrary to the assumption f75|) . Note that for a given network connectivity, the condition f92|) is satisfied for any 
Z > l c where l c is the diameter of the underlying network, the longest directed connection path between any two 
oscillators in the network. The diameter is maximal, l c — N — 1, assumed for a ring of N oscillators. 
Hence, after at most l c periods, the norm of the perturbation vector decreases, 

\\6(IT)\\ < \\S\\ (95) 

for I > l c , Since ||<5(mZ c T)|| is strictly monotonically decreasing with m 6 N and bounded below by zero, the limit 

S°° := lim 8{ml c T) (96) 

m — >oc 

exists. If S°° ciVi the vector norm would be reduced as implied by the above considerations. Thus we find that 

S°° = ciVi (97) 

because the uniform components civi of the original perturbation S does not change under the dynamics (see Eq. i[73|) ). 
This completes the demonstration that the synchronous state is asymptotically stable for any strongly connected 
network of inhibitorily coupled oscillatory units. 



VIII. SUMMARY AND DISCUSSION 



We analyzed the stability of synchronous states in arbitrarily connected networks of pulse-coupled oscillators. For 
generally structured networks, the intricate problem of multiple stability operators arises. We analyzed this multi- 
operator problem exactly. For both inhibitory and excitatory couplings, we determined analytical bounds on the 
eigenvalues of all operators. Given the multi-operator property, stability cannot be deduced by considering the 
eigenvalues only. We therefore completed the stability analysis on two levels For networks with inhibitory couplings 
(eij < 0) we found plain (Lyapunov) stability of the synchronous state; under mild constraints (strong connectivity of 
the network), graph theoretical arguments show that it is also asymptotically stable, independent of the parameters 
and the details of the connectivity structure. 

Let us point out that the stability results obtained here are valid for arbitrarily large finite networks of general 
connectivity. In contrast to ordinary stability problems, the eigenvalues of the first order operators do not directly 
determine the stability here, because different linear operators may act subsequently on a given perturbation vector. 
Thus the eigenvalues and eigenvectors of an individual matrix alone do not define the local dynamics of the system. 
For network dynamical systems of pulse-coupled units that exhibit more complex cluster states, with two or more 
groups of synchronized units, an analogous multi-operator problem arises [13, Hil Ell] ■ Methods to bound their spectra 
and the graph theoretical approach presented above to prove asymptotic stability are applicable to these cluster states 
in a similar way, provided the stability problem is originally formulated in an event-based manner. Further studies 
(e.g., [29]) show that the stable synchrony in complex networks is not restricted to the class of models considered 
here. Together with other works, e.g. [Ill, [H, HH, this also indicates robustness of the synchronous state against 
structural perturbations, i.e. conditions necessary to ensure that the local dynamics stays qualitatively similar if the 
time evolution of the model is weakly modified. 
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It has been hypothesized (see e.g. [-481 . [501) that the experimentally observed precisely timed firing patterns in 
otherwise irregular neural activity might be generated dynamically by the cerebral cortex. As a by-product, our 
results clearly demonstrate that states in which units fire in a temporally highly regular fashion and states with 
irregular asynchronous activity may be coexisting attractors of the same recurrent network, cf. [38]. This already 
applies to random networks. More specifically structured networks may possess a large variety of dynamical states in 
which firing times are precisely coordinated. A promising direction for future research is thus to adopt the methods 
developed here for investigating the stability of such states. In particular, methods adapted from the ones developed 
here may reveal dynamical features of networks in which the coupling strengths are not only structured but highly 
heterogeneous and precise temporal firing patterns occur in place of a simple synchronous state [29l I34L l36| . More 
generally, our methods may help to understand properties of stability and robustness in various network dynamical 
systems with pulsed interactions. They may thus be relevant not only for networks of nerve cells in the brain coupled 
via chemical synapses but also for cells in heart tissue coupled electrically and even for populations of fireflies and 
chirping crickets that interact by sending and receiving light pulses and sound signals, respectively. 

Acknoledgements: We thank Theo Geisel, Michael Denker and Raoul-Martin Memmesheimer for fruitful discussions. 
This work has been partially supported by the Ministry of Education and Research (BMBF) Germany under Grant 
No. 01GQ07113 and, via the Bernstein Center for Computational Neuroscience (BCCN) Gottingen, under Grant No. 
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IX. APPENDIX: EXACT DERIVATION OP THE EXPANSION ((34)) 

We prove a generalization of the first order expansion |34]) used above, 

m 

0i,m = Oi,m + /] Pi,n-l,m A,n for m £ {1, . . . , fcj, 



(98) 



by induction over m for arbitrary fixed i £ {1, . . . , N}. As in the main text, here the statement x = y means that 
x = y + J^i n O(Di.n) as a ^ A,n — * 0- The quantities /?j, m were defined in Table [J . The quantities 1(98]) appear in 
the derivation of the map lf32]) for m = h. As an extension of the notation in section llVl we denote here 



= U- l \U{T) + Y j e ij \ for me{0,...,ki}, 



n=l 



(99) 



and 



A,n = A i>n _i - A ijn for ne{l,...,h}, 
U'(oi, n ) 



Pi 



for < n < m < ki 



(100) 



(101) 



Here the prime denotes the derivative of the potential function U with respect to its argument. The latter definition 
implies the identity 



Pi,n,l Pi,l,m — Pi,n,m 



For the case m = 1, the induction basis, expression (|98| holds because 

= U-\U(T + D hl )+e l3l ) 
= U-\U{r)+e lh ) 



' d_ 

dD 

Ui.i + 

CtiA + 



U-\U{T + D) + e in ) 
U'(r) 

U'(U^(U(T)+e in )) 

tf'Ko) 



D 



i.l 



D=0 



U'{at,i) 

Oii,l + PifiA A,l 



A,i 



(102) 



(103) 
(104) 



(105) 

(106) 
(107) 
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where the third last equality holds because the derivatives of inverse functions are given by d/dxU (x) 
l/U'(U^ 1 (x)). The induction step, m i— ► m+ 1, is proven by 



A,m+1 — U 1 {U{(3i t m + Di^ n+ i) + £ij m + 1 ) 

= u- 1 (u ( 



= U^ 1 (U(ai, m ) +£ij m+1 ) 
d 



dD 



U- 1 {U(a hm + D)+£ lJm+1 ) 



D=0 \„ = 



^ ^ Pi.n—\,ni ^i,n -{- D.[ m + i 



n=l 



U 1 (U (a iym ) + £ij m+1 ) + TT/( ^ ""^ ( ^2pi,n-l,m D i>n + A,m+1 



f '(ai,m+l) 



\n=l 



"bP^rr^m+l I ^ ^ Pi,n — l,m Dj,n ~b 

\n=l / 

m+l 

~t~ ^ ^ Pi,n — l,m+l ^i.n 
n=l 



(108) 
(109) 
(110) 

(111) 
(112) 
(113) 



where we used J7(ai jTO ) + £ij m+1 = f7(a^ m+ i) in the second last step and the identity (|102|) in the last step. This 
completes the proof of the first order expansion lj98j) . 



Note that because of the normalization Yln=i E ijn ~ e *i = e f° r a ^ *i ^ ne Q uan tity &i,ki = U 1 {U(t)+e) = a 

is independent of i. In addition, Pi, n ,ki = Pi,™ for all i and all n. Hence, theorem (|98|) in the case m = ki yields 
expression (|34|l . 
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